/* External definitions for single-server queueing system. */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>  /* Header file for random-number generator. */

#define Q_LIMIT 100  /* Limit on queue length. */
#define BUSY      1  /* Mnemonics for server's being busy */
#define IDLE      0  /* and idle. */
#define a (double) 16807.0
#define m (double) 2147483647.0
#define q (double) 127773.0
#define r (float) 2836.0

int   next_event_type, num_custs_delayed, num_delays_required, num_events,
num_in_q, server_status, job_gt_45, count;
float area_num_in_q, area_server_status, mean_interarrival, mean_service,
sim_time, time_arrival[Q_LIMIT + 1], time_last_event, time_next_event[3],
total_of_delays, my_random, time_required, service_times[Q_LIMIT+1], time_q_gt1, avg_time_in_q[10], avg_in_sys[10], gt_45[10];
FILE  *infile, *outfile;
double hi, test, lo, *seed1, *seed2;
char service, type;

void  initialize(void);
void  timing(void);
void  arrive(void);
void  depart(void);
void  report(void);
void  update_time_avg_stats(void);
float expon(float mean, double*);
float randomvar(double*);


main()  /* Main function. */
{
	/* Specify the number of events for the timing function. */
	double seed1array[10] = {1806794933,1,2,3,4,5,6,7,8,9};
	double seed2array[10] = {1155192169,11111,22222,33333,44444,55555,66666,77777,88888,99999};
	outfile = fopen("mm1.out", "w");
	num_events = 2;

	/* Specify the starting Seed */
	char terminating_condition;
	do
	{
		printf("Enter if service times are M or D\n");
		scanf("%c", &service);
		if(service != 'M' && service != 'D')
		{
			printf("Incorrect input, try again\n");
		}
	}
	while(service != 'M' && service != 'D');
	if(service == 'M'){
		do
		{
			printf("Enter (F)CFS or (S)PT\n");
			scanf("\n%c", &type);
			if(type != 'F' && type != 'S')
			{
				printf("Incorrect input, try again\n");
			}
		}
		while(type != 'F' && type != 'S');
	}

	int done = 0;
	seed1 = malloc(sizeof(double));
	seed2 = malloc(sizeof(double));
	int i;

	do
	{
		printf("Enter terminating condition time(t) or people(p)\n");
		scanf("\n%c", &terminating_condition);
		if(terminating_condition == 't')
		{
			printf("Enter terminating time in min\n");
			scanf("%f", &time_required);
			done = 1;
		}
		else if(terminating_condition == 'p')
		{
			printf("Enter terminating number of people\n");
			scanf("%d", &num_delays_required);
			time_required = -1;
			done = 1;
		}
		else
		{
			printf("Incorrect format try again");
		}
	}
	while(done != 1);
	/* Read input parameters. */
	fprintf(outfile, "Average delay in queue		Average number in system		Fraction of Jobs that spent more than 4.5 minutes in the system\n\n");
	mean_interarrival = 5;
	mean_service = 4;
	/* Write report heading and input parameters. */
	for(count = 0; count < 10; count++)
	{
		*seed1 = seed1array[count];
		*seed2 = seed2array[count];
		/* Initialize the simulation. */
		initialize();

		/* Run the simulation while more delays are still needed. */
		if(terminating_condition == 'p')
		{
			while (1) {
				/* Determine the next event. */
				timing();
				/* Update time-average statistical accumulators. */

				update_time_avg_stats();
				/* Invoke the appropriate event function. */
				if(num_custs_delayed < num_delays_required || next_event_type == 1)
				{
					switch (next_event_type) {
					case 1:
						arrive();
						break;
					case 2:
						depart();
						break;
					}
				}
				else
				{
					break;
				}
			}
		}
		else
		{
			while (1) {
				/* Determine the next event. */

				timing();
				if(sim_time > time_required)
				{
					sim_time = time_required;
				}
				update_time_avg_stats();
				/* Update time-average statistical accumulators. */
				if(sim_time < time_required) {


					/* Invoke the appropriate event function. */

					switch (next_event_type) {
					case 1:
						arrive();
						break;
					case 2:
						depart();
						break;
					}
				}
				else
				{
					num_custs_delayed--;
					break;
				}
			}
		}
		report();
	}

	/* Invoke the report generator and end the simulation. */



	fclose(outfile);

	return 0;
}


void initialize(void)  /* Initialization function. */
{
	/* Initialize the simulation clock. */
	int i;
	sim_time = 0.0;

	/* Initialize the state variables. */

	server_status   = IDLE;
	num_in_q        = 0;
	time_last_event = 0.0;

	/* Initialize the statistical counters. */
	job_gt_45          = 0;
	num_custs_delayed  = 0;
	total_of_delays    = 0.0;
	area_num_in_q      = 0.0;
	area_server_status = 0.0;
	time_q_gt1		   = 0.0;
	/* Initialize event list.  Since no customers are present, the departure
       (service completion) event is eliminated from consideration. */

	time_next_event[1] = sim_time + expon(mean_interarrival, seed1);
	time_next_event[2] = 1.0e+30;
	for(i = 0; i < Q_LIMIT; i++)
	{
		service_times[i] = 0;
	}
}


void timing(void)  /* Timing function. */
{
	int   i;
	float min_time_next_event = 1.0e+29;

	next_event_type = 0;

	/* Determine the event type of the next event to occur. */

	for (i = 1; i <= num_events; ++i)
		if (time_next_event[i] < min_time_next_event) {
			min_time_next_event = time_next_event[i];
			next_event_type     = i;
		}

	/* Check to see whether the event list is empty. */

	if (next_event_type == 0) {

		/* The event list is empty, so stop the simulation. */

		fprintf(outfile, "\nEvent list empty at time %f", sim_time);
		exit(1);
	}

	/* The event list is not empty, so advance the simulation clock. */
	sim_time = min_time_next_event;
}


void arrive(void)  /* Arrival event function. */
{
	float delay;

	/* Schedule next arrival. */
	time_next_event[1] = sim_time + expon(mean_interarrival, seed1);
	/* Check to see whether server is busy. */

	if (server_status == BUSY) {

		/* Server is busy, so increment number of customers in queue. */

		++num_in_q;
		/* Check to see whether an overflow condition exists. */

		if (num_in_q > Q_LIMIT) {

			/* The queue has overflowed, so stop the simulation. */

			fprintf(outfile, "\nOverflow of the array time_arrival at");
			fprintf(outfile, " time %f", sim_time);
			exit(2);
		}

		if(type == 'S' && service == 'M')
		{
			float new_service_time =  expon(mean_service, seed2);
			int i;
			if(service_times[1] == 0)
			{
				service_times[1] = new_service_time;
				time_arrival[1] = sim_time;
			}
			else
			{
				for(i = Q_LIMIT; i > 1; i--)
				{
					service_times[i] = service_times[i-1];
					time_arrival[i] = time_arrival[i-1];
					if(service_times[i] < new_service_time && service_times[i] != 0)
					{
						service_times[i] = new_service_time;
						time_arrival[i] = sim_time;
						break;
					}
					else if(i == 2)
					{
						time_arrival[1] = sim_time;
						service_times[1] = new_service_time;
					}
				}
			}
		}
		else
		{
			time_arrival[num_in_q] = sim_time;
		}


		/* There is still room in the queue, so store the time of arrival of the
           arriving customer at the (new) end of time_arrival. */


	}

	else {

		/* Server is idle, so arriving customer has a delay of zero.  (The
           following two statements are for program clarity and do not affect
           the results of the simulation.) */

		delay            = 0.0;
		total_of_delays += delay;

		/* Increment the number of customers delayed, and make server busy. */

		++num_custs_delayed;
		server_status = BUSY;

		/* Schedule a departure (service completion). */

		if(service == 'M')
		{
			time_next_event[2] = sim_time + expon(mean_service, seed2);

			//printf("Customer Delayed no delay: %d, finished: %f\n", num_custs_delayed, time_next_event[2]);
		}
		else if(service == 'D')
		{
			time_next_event[2] = sim_time + mean_service;
		}
		if(time_next_event[2] - sim_time + delay > 4.5)
		{
			if(time_required > time_next_event[2])
			{
				job_gt_45++;
			}
			else if(time_required < 0)
			{
				job_gt_45++;
			}
		}
	}
}


void depart(void)  /* Departure event function. */
{
	int   i;
	float delay;

	/* Check to see whether the queue is empty. */

	if (num_in_q == 0) {

		/* The queue is empty so make the server idle and eliminate the
           departure (service completion) event from consideration. */

		server_status      = IDLE;
		time_next_event[2] = 1.0e+30;
	}

	else {

		/* The queue is nonempty, so decrement the number of customers in
           queue. */

		--num_in_q;

		/* Compute the delay of the customer who is beginning service and update
           the total delay accumulator. */

		delay            = sim_time - time_arrival[1];
		total_of_delays += delay;

		/* Increment the number of customers delayed, and schedule departure. */

		++num_custs_delayed;
		if(service == 'M')
		{
			if(type == 'S')
			{
				time_next_event[2] = sim_time + service_times[1];
				//printf("Customer Delayed: %d, finished: %f\n", num_custs_delayed, time_next_event[2]);
				for(i = 1; i < Q_LIMIT; i++)
				{
					service_times[i] = service_times[i+1];
				}
				service_times[Q_LIMIT] = 0;
			}
			else
			  {

			    time_next_event[2] = sim_time + expon(mean_service, seed2);
			    //printf("Customer Delayed: %d, finished: %f\n", num_custs_delayed, time_next_event[2]);
			  }
		}
		else if(service == 'D')
		{
			time_next_event[2] = sim_time + mean_service;
		}
		if(time_next_event[2] - sim_time + delay > 4.5)
		{
			if(time_required > time_next_event[2])
			{
				job_gt_45++;
			}
			else if(time_required < 0)
			{
				job_gt_45++;
			}
		}
		/* Move each customer in queue (if any) up one place. */

		for (i = 1; i <= num_in_q; ++i)
			time_arrival[i] = time_arrival[i + 1];
	}
}


void report(void)  /* Report generator function. */
{
	/* Compute and write estimates of desired measures of performance. */
	float server_util = area_server_status/sim_time;
	float avg_in_system = (area_server_status+area_num_in_q)/sim_time;
	fprintf(outfile, "%11.3f%30.3f%22.3f\n", total_of_delays / (float)num_custs_delayed, avg_in_system, (float)job_gt_45/(float)num_custs_delayed);
	avg_in_sys[count] = avg_in_system;
	avg_time_in_q[count] = total_of_delays /(float)num_custs_delayed;
	gt_45[count] = (float)job_gt_45/(float)num_custs_delayed;
	if(count == 9)
	{
		int i;
		float mean_in_sys, mean_time_q, mean_gt, var_sys, var_q, var_gt;
		mean_in_sys = 0;
		mean_time_q = 0;
		mean_gt = 0;
		var_sys = 0;
		var_q = 0;
		var_gt = 0;
		for(i = 0; i < 10; i++)
		{
			mean_in_sys += avg_in_sys[i];
			mean_time_q += avg_time_in_q[i];
			mean_gt += gt_45[i];
		}
		mean_in_sys = mean_in_sys/10;
		mean_time_q = mean_time_q/10;
		mean_gt = mean_gt/10;
		fprintf(outfile, "Xbar(n)%6.3f%30.3f%22.3f\n", mean_time_q, mean_in_sys, mean_gt);
		for(i = 0; i < 10; i++)
		{
			var_sys += pow(avg_in_sys[i]-mean_in_sys,2.0)/9.0;
			var_q += pow(avg_time_in_q[i]-mean_time_q,2.0)/9.0;
			var_gt += pow(gt_45[i]-mean_gt,2.0)/9.0;
		}
		fprintf(outfile, "S^2(n)%6.3f%30.3f%22.3f\n", var_q, var_sys, var_gt);
		fprintf(outfile, "sqrt(S^2(n)/n)%3.3f%30.3f%22.3f\n", sqrt(var_q/10), sqrt(var_sys/10), sqrt(var_gt/10));

	}
}


void update_time_avg_stats(void)  /* Update area accumulators for time-average
                                     statistics. */
{
	float time_since_last_event;

	/* Compute time since last event, and update last-event-time marker. */

	time_since_last_event = sim_time - time_last_event;
	time_last_event       = sim_time;

	/* Update area under number-in-queue function. */

	area_num_in_q      += num_in_q * time_since_last_event;
	if(num_in_q > 1)
	{
		time_q_gt1 += time_since_last_event;
	}
	/* Update area under server-busy indicator function. */

	area_server_status += server_status * time_since_last_event;
}


float expon(float mean, double *ex_seed)  /* Exponential variate generation function. */
{
	/* Return an exponential random variate with mean "mean". */

	return -mean * log(randomvar(ex_seed));
}

float randomvar(double *rand_seed)
{      
	hi=floor(*rand_seed/q);
	lo=*rand_seed-q*hi;
	test=a*lo-r*hi;
	if (test>0.0)
		*rand_seed=test;
	else
		*rand_seed=test+m;
	my_random=*rand_seed/m;
	return my_random;
}
